##############################
# Storage blog - March, 2018 #
##############################

rm(list=ls())
root <- 'C:/Users/steve/OneDrive/Blogs/storage mar 2018/' # Set this to the local file-path where the EIA data is stored
library(readxl)
library(dplyr)
library(ggplot2)

# Load storage data from EIA860
a <- read_excel(paste0(root, '3_4_Energy_Storage_Y2016.xlsx'), skip=1)

# Make application dummies
for(i in 1:nrow(a)){
  for(j in 27:37){
    if(a[i,j]=='Y' & !is.na(a[i,j])){a[i,j] <- 1} 
    else{a[i,j] <- 0}
  }
}

# Merge location data from EIA860
b <- read_excel(paste0(root, '2___Plant_Y2016.xlsx'), skip=1)
b <- b[,c('Plant Code', 'Street Address', 'City','Zip', 'Latitude', 'Longitude')]
a <- merge(a, b, by='Plant Code')
a <- subset(a, !a$State%in%c('AK','HI')) 
a <- a[c(1:13, 15:nrow(a)),] 
a[13, 'Nameplate Capacity (MW)'] <- 280 # Combine twos solar thermal generators on same site into one plant

# Make map
a <- arrange(a, desc(`Nameplate Capacity (MW)`))
usa <- map_data('state')
ggplot() + geom_polygon(data=usa, aes(x=long, y=lat, group=group), fill='ivory', color='grey70') + coord_fixed(1.4) + theme_void() + 
  geom_point(data=subset(a,!a$State%in%c('AK','HI')), aes(x=Longitude, y=Latitude, size=`Nameplate Capacity (MW)`, fill=Technology), shape=21, color='black') +
  scale_size(limits=c(0.5,280)) + scale_fill_manual(values=c('salmon','steelblue','darkolivegreen3','darkorchid2')) +
  guides(fill = guide_legend(override.aes = list(size=4))) +
  theme(legend.title = element_text(size=10), legend.text = element_text(size=10)) + labs(caption = "Energy storage deployment, 2016.") +
  theme(plot.caption = element_text(size=14))
ggsave(paste0(root, 'map.png'), width=9, height=6)

# Plot annual growth over time
a <- a[which(!a$Technology%in%c('Natural Gas with Compressed Air Storage','Solar Thermal with Energy Storage')),] # Remove these two technologies from chart
d <- data.frame(Year=numeric(0), `Nameplate Capacity (MW)`=numeric(0))
for(yr in 2008:2016){
  e <- data.frame(Year=numeric(1))
  e$Year <- yr
  e$`Nameplate Capacity (MW)` <- sum(a[which(a$`Operating Year`<=yr),'Nameplate Capacity (MW)'])
  d <- rbind(d,e)
}
d$Year <- factor(d$Year)
ggplot(data=d, aes(x=Year, y=`Nameplate Capacity (MW)`)) + geom_bar(stat='identity', fill='steelblue') + 
  theme_bw() + theme(text=element_text(size=18)) 
ggsave(paste0(root, 'growth.png'), width=6, height=4)


# Make services chart
c <- data.frame(Service=character(11), stringsAsFactors = F)
c[1,'Service'] <- 'Arbitrage'
c[1, 'Nameplate Capacity (MW)'] <- sum(a[which(a$Arbitrage==1),'Nameplate Capacity (MW)'])
c[2,'Service'] <- 'Frequency Regulation'
c[2, 'Nameplate Capacity (MW)'] <- sum(a[which(a$`Frequency Regulation`==1),'Nameplate Capacity (MW)'])
c[3,'Service'] <- 'Load Following'
c[3, 'Nameplate Capacity (MW)'] <- sum(a[which(a$`Load Following`==1),'Nameplate Capacity (MW)'])
c[4,'Service'] <- 'Ramping / Spinning Reserve'
c[4, 'Nameplate Capacity (MW)'] <- sum(a[which(a$`Ramping / Spinning Reserve`==1),'Nameplate Capacity (MW)'])
c[5,'Service'] <- 'Co-Located Renewable Firming'
c[5, 'Nameplate Capacity (MW)'] <- sum(a[which(a$`Co-Located Renewable Firming`==1),'Nameplate Capacity (MW)'])
c[6,'Service'] <- 'Transmission and Distribution Deferral'
c[6, 'Nameplate Capacity (MW)'] <- sum(a[which(a$`Transmission and Distribution Deferral`==1),'Nameplate Capacity (MW)'])
c[7,'Service'] <- 'System Peak Shaving'
c[7, 'Nameplate Capacity (MW)'] <- sum(a[which(a$`System Peak Shaving`==1),'Nameplate Capacity (MW)'])
c[8,'Service'] <- 'Load Management'
c[8, 'Nameplate Capacity (MW)'] <- sum(a[which(a$`Load Management`==1),'Nameplate Capacity (MW)'])
c[9,'Service'] <- 'Voltage or Reactive Power Support'
c[9, 'Nameplate Capacity (MW)'] <- sum(a[which(a$`Voltage or Reactive Power Support`==1),'Nameplate Capacity (MW)'])
c[10,'Service'] <- 'Backup Power'
c[10, 'Nameplate Capacity (MW)'] <- sum(a[which(a$`Backup Power`==1),'Nameplate Capacity (MW)'])
c[11,'Service'] <- 'Excess Wind and Solar Generation'
c[11, 'Nameplate Capacity (MW)'] <- sum(a[which(a$`Excess Wind and Solar Generation`==1),'Nameplate Capacity (MW)'])
c <- arrange(c, desc(`Nameplate Capacity (MW)`))
c$Service <- factor(c$Service, levels=rev(c('Frequency Regulation','Ramping / Spinning Reserve','Load Following', 'Voltage or Reactive Power Support',
                                            'System Peak Shaving','Load Management','Excess Wind and Solar Generation', 'Arbitrage',   'Transmission and Distribution Deferral',
                                           'Backup Power','Co-Located Renewable Firming')))

ggplot(data=c, aes(x=Service, y=`Nameplate Capacity (MW)`)) + geom_bar(stat='identity', fill='steelblue') + coord_flip() + theme_bw()+ theme(text=element_text(size=12))
ggsave(paste0(root, 'services.png'), width=6, height=4)

# Analyze active storage
gen <- read_excel(paste0(root, 'EIA923_Schedules_2_3_4_5_M_12_2016_Final_Revision.xlsx'), skip=5)
names(gen)[14] <- 'prime_mover'
gen <- subset(gen, gen$prime_mover %in% c('BA','FW'))
names(gen)[96] <- 'netgen'
gen <- gen[, c('Plant Id', 'netgen')]
gen <- rename(gen, `Plant Code` = `Plant Id`)
a <- merge(a,gen, by='Plant Code')
a[which(a$netgen==0),'Technology'] <- 'Inactive'
( activemw <- sum(a[which(a$netgen!=0),'Nameplate Capacity (MW)']) )
( inactivemw <- sum(a[which(a$netgen==0),'Nameplate Capacity (MW)']) )
( activeshare <- activemw / (activemw + inactivemw) )

ggplot() + geom_polygon(data=usa, aes(x=long, y=lat, group=group), fill='ivory', color='grey70') + coord_fixed(1.4) + theme_void() + 
  geom_point(data=subset(a,!a$State%in%c('AK','HI')), aes(x=Longitude, y=Latitude, size=`Nameplate Capacity (MW)`, fill=Technology), shape=21, color='black') +
  scale_size(limits=c(0.5,50)) + scale_fill_manual(values=c('salmon','steelblue','black')) +
  guides(fill = guide_legend(override.aes = list(size=4))) +
  theme(legend.title = element_text(size=10), legend.text = element_text(size=10)) + labs(caption = "Active vs. inactive storage facilities, 2016.") +
  theme(plot.caption = element_text(size=14))
ggsave(paste0(root, 'activemap.png'), width=9, height=6)
